﻿using System;
using System.Collections;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Drawing.Imaging;
using System.Drawing;
using System.Windows.Forms;

namespace SARD.Classes
{
    public class dicomimage
    {
        DicomDecoder dd;
        byte[] pixels8;
        ushort[] pixels16;
        int Width;
        int Height;
        int bitDepth;
        int max, min;
        Bitmap bmp;
        int odlM, odlm;
        bool zeroiswhite;

        /// <summary>
        /// Basic class constructor. It will set the pixel list to null (8 and 16 bits)
        /// </summary>
        public dicomimage() 
        {
            pixels16 = null;
            pixels8 = null;
            zeroiswhite = false;
        }

        /// <summary>
        /// Basic class constructor. It will create a 16 bits pixel list
        /// </summary>
        public dicomimage(ushort[] dicominfo)
        {
            pixels16 = new ushort[dicominfo.Length];
            pixels16 = dicominfo;
            pixels8 = null;
            zeroiswhite = false;
        }

        /// <summary>
        /// Basic class constructor. It will create a 8 bits pixel list
        /// </summary>
        public dicomimage(byte[] dicominfo)
        {
            pixels8 = new byte[dicominfo.Length];
            pixels8 = dicominfo;
            pixels16 = null;
            zeroiswhite = false;
        }

        /// <summary>
        /// Class dicom image constructor.
        /// It's the point where the image decoding starts.
        /// It will call the proper functions to read and display the image.
        /// </summary>
        /// <param name="file"></param>
        public dicomimage(string file)
        {
            pixels16 = null;
            pixels8 = null;
            dd = new DicomDecoder();
            zeroiswhite = false;
            ReadAndDisplayDicomFile(file);

        }


        /// <summary>
        /// Return the depth of the image (dicom).
        /// </summary>
        /// <returns>
        /// Depth as a int
        /// </returns>
        public int bitdepth() 
        {
            return this.bitDepth;
        }

        /// <summary>
        /// Returns the last used value of minimum (for displaying propouses).
        /// </summary>
        /// <returns>
        /// Old minimum as a int
        /// </returns>
        public int getoldm()
        {
            return this.odlm;

        }

        /// <summary>
        /// Returns the last used value of maximum (for displaying propouses).
        /// </summary>
        /// <returns>
        /// Old maximum as a int
        /// </returns>
        public int getoldM()
        {
            return this.odlM;

        }

        /// <summary>
        /// Returns the value of minimum (for displaying propouses).
        /// </summary>
        /// <returns>
        /// Minimum as a int
        /// </returns>
        public int getmin()
        {
            return this.min;
        }

        /// <summary>
        /// Returns the value of maximum (for displaying propouses).
        /// </summary>
        /// <returns>
        /// Maximum as a int
        /// </returns>
        public int getmax()
        {
            return this.max;
        }


        /// <summary>
        /// Returns the value of the image width
        /// </summary>
        /// <returns>
        /// Width as a int
        /// </returns>
        public int getwidth()
        {
            return this.Width;
        }

        /// <summary>
        /// Returns the value of the image height
        /// </summary>
        /// <returns>
        /// Height as a int
        /// </returns>
        public int getheight()
        {
            return this.Height;
        }

        public void changereadcolor(bool type)
        {
            zeroiswhite = type;
        }

        /// <summary>
        /// Returns a bitmap image with one is the coverted dicom image.
        /// This is a high level function that will find the correct parameters in use and proced to the best answer
        /// </summary>
        /// <param name="m">
        /// It's the smallest value of the desired value range to be displayed.
        /// </param>
        /// <param name="M">
        /// It's the biggest value of the desired value range to be displayed.
        /// </param>
        /// <returns></returns>
        public Bitmap display(int m, int M) 
        {
            this.odlm = m;
            this.odlM = M;

            if (bitDepth == 8)
            {
                changescaleofpixelsdisplayed8(m, M);
                return bmp;
            }
            else
            {
                changescaleofpixelsdisplayed16(m, M);
                return bmp;
            }

        }

        /// <summary>
        /// Returns a bitmap image with one is the coverted dicom image.
        /// This is a high level function that will find the correct parameters in use and proced to the best answer
        /// </summary>
        /// <returns></returns>
        public Bitmap display()
        {
            if (bitDepth == 8)
            {
                setmaxandmin8();
                changescaleofpixelsdisplayed8();
                return bmp;
            }
            else
            {
                setmaxandmin16();
                changescaleofpixelsdisplayed16();
                return bmp;
            }
        }

        /// <summary>
        /// Returns the tags of the associated dicom file
        /// </summary>
        /// <returns></returns>
        public List<string> tags()
        {
            return this.dd.dicomInfo;
        }

        /// <summary>
        /// It acts as a image filter.
        /// </summary>
        /// <returns></returns>
        public Bitmap invertimage()
        {
            if (bitDepth == 8)
            {
                invertpixels8();
                changescaleofpixelsdisplayed8();
            }
            else
            {
                invertpixels16();
                changescaleofpixelsdisplayed16();
            }

            return bmp;

        }

        /// <summary>
        /// Calculate the average of the values of a desired rectangle inside the image.
        /// </summary>
        /// <param name="calc">
        /// Coordinates of the desired piece of image where the funcion will work.
        /// </param>
        /// <returns></returns>
        public double average(Rectangle calc) 
        {
            if (bitDepth == 8)
                return averagevalue8(calc);

            else
                return averagevalue16(calc);
        }

        /// <summary>
        /// Calculate the stddev of the values of a desired rectangle inside the image.
        /// </summary>
        /// <param name="calc">
        /// Coordinates of the desired piece of image where the funcion will work.
        /// </param>
        /// <returns></returns>
        public double stddev(Rectangle calc)
        {
            if (bitDepth == 8)
                return stddev8(calc,  averagevalue8(calc));

            else
                return stddev16(calc, averagevalue16(calc));

        }


        /// <summary>
        /// It is desing to send the min and the max value of the range to the user
        /// </summary>
        /// <param name="calc">
        /// Coordinates of the desired piece of image where the funcion will work.
        /// </param>
        /// <returns></returns>
        public Point maxandmin(Rectangle calc) 
        {
            if (bitDepth == 8)
               return maxandmin8(calc);
            else
                return maxandmin16(calc);
        }
        

        /*
        public Point maxandminoveraline(Point p1, Point p2)
        {
            if (bitDepth == 8)
                return maxandmin8overaline(p1, p2);
            else
                return maxandmin16overaline(p1, p2);
        }
        */


        /// <summary>
        /// Apply the proper functions to star the unifformidity calculous.
        /// </summary>
        /// <param name="size">
        /// It is the desired size of the square size (mimification image).
        /// </param>
        public void uniformidity(int size) 
        {
            if (bitDepth == 8)
            {
                acumulateimage8(size);
                filter16(size);
                removebadpixels(size);

            }

            else 
            {
                acumulateimage16(size);
                filter16(size);
                removebadpixels(size);

            }

        }

        /// <summary>
        /// Returns the integral unifformidy calculated over a rectangle surface.
        /// </summary>
        /// <param name="calc">
        /// Coordinates of the desired piece of image where the funcion will work.
        /// </param>
        /// <returns></returns>
        public double IntegralUnifform(Rectangle calc) 
        {
            Point calculate = maxandmin(calc);
            return (double)(calculate.Y - calculate.X) / (calculate.Y + calculate.X);
        }

        /// <summary>
        /// Returns the differential unifformidy calculated over a rectangle surface.
        /// </summary>
        /// <param name="calc">
        /// Coordinates of the desired piece of image where the funcion will work.
        /// </param>
        /// <returns></returns>
        public double DifferentialUnifform(Rectangle calc)
        {
            return unifordiff16(calc);
        }

        /// <summary>
        /// Returns the value of a specific position (for displaying propouses).
        /// </summary>
        /// <param name="i">
        /// The line position of the value.
        /// </param>
        /// <param name="j">
        /// The columm position of the value.
        /// </param>
        /// <returns></returns>
        public int returnvalue(int i, int j) 
        {
            if (i >= 0 && i <= this.Width && j >= 0 && j <= this.Height)
            {
                if (bitDepth == 8)
                    return pixels8[i + j * Width];

                if (bitDepth == 16)
                    return pixels16[i + j * Width];

                else
                    return 0;
            }
            else
                return 0;

        }

        /// <summary>
        /// Populates the histogram created by a image piece.
        /// </summary>
        /// <param name="calc">
        /// It's the area where the calculation will take place.
        /// </param>
        /// <returns>
        /// Returns an array with the histogram data.
        /// </returns>
        public double [] generatehistogram(Rectangle calc) 
        {
            if (calc.Width > 0)
            {
                Point p = this.maxandmin(calc);

                int lenght = p.Y - p.X + 1;
                double[] histogramarray = new double[lenght];

                for (int i = 0; i < lenght; i++)
                    histogramarray[i] = 0;


                if (bitDepth == 8)
                    for (int j = 0; j < calc.Height; j++)
                        for (int i = calc.X; i <= calc.X + calc.Width; i++)
                            histogramarray[pixels8[i + (j + calc.Y) * Width] - p.X]++;

                if (bitDepth == 16)
                    for (int j = 0; j < calc.Height; j++)
                        for (int i = calc.X; i <= calc.X + calc.Width; i++)
                            histogramarray[pixels16[i + (j + calc.Y) * Width] - p.X]++;


                return histogramarray;
            }
            else 
                return null;
        }


        /// <summary>
        /// Reads the dicom image and create an array of data. It also sets the image parameters.
        /// </summary>
        /// <param name="fileName"></param>
        private void ReadAndDisplayDicomFile(string fileName)
        {
            dd.DicomFileName = fileName;
            bool result = dd.dicomFileReadSuccess;
            if (result == true)
            {
                Width = dd.width;
                Height = dd.height;
                bitDepth = dd.bitsAllocated;


                if (bitDepth == 8)
                {
                    pixels8 = null;

                    //The bitmap image that has been read will apear here.
                    dd.GetPixels8(ref pixels8);

                    setmaxandmin8();


                    // Perform window level using the image min and max values;
                    //I will also create another list to ensure that the base values will not be changed.
                    changescaleofpixelsdisplayed8( min, max);

                }

                if (bitDepth == 16)
                {
                    pixels16 = null;


                    dd.GetPixels16(ref pixels16);

                    setmaxandmin16();

                    changescaleofpixelsdisplayed16(min, max);

                }
            }

            else
            {
                if (dd.dicmFound == false)
                {
                    MessageBox.Show("Este arquivo não é válido.");
                }
                else if (dd.dicomDir == true)
                {
                    MessageBox.Show("Este arquivo parece com um arquivo DICOM, mas não contém uma imagem.");
                }
                else
                {
                    MessageBox.Show("Não foi possível ler o arquivo solicitado\n" +
                        "No entanto será possível visualizar as Tasgs do arquivo.");
                }

                // Show a plain grayscale image

                bmp = new Bitmap(1024, 1024, PixelFormat.Format24bppRgb);

                BitmapData bmd = bmp.LockBits(new Rectangle(0, 0, 1024, 1024),
                    System.Drawing.Imaging.ImageLockMode.ReadOnly, bmp.PixelFormat);

                unsafe
                {
                    int pixelSize = 3;
                    int i, j, j1;
                    byte bite;

                    for (i = 0; i < bmd.Height; ++i)
                    {
                        byte* row = (byte*)bmd.Scan0 + (i * bmd.Stride);

                        for (j = 0; j < bmd.Width; ++j)
                        {

                            bite = 224;
                            j1 = j * pixelSize;
                            row[j1] = bite;            // Red
                            row[j1 + 1] = bite;        // Green
                            row[j1 + 2] = bite;        // Blue
                        }
                    }
                } // end unsafe      
                bmp.UnlockBits(bmd);

            }

        }


        /// <summary>
        /// Read the image pixels and create a bitmap for displaying the image.
        /// </summary>
        /// <param name="m">
        /// It's the smallest value of the desired value range to be displayed.
        /// </param>
        /// <param name="M">
        /// It's the biggest value of the desired value range to be displayed.
        /// </param>
        private void changescaleofpixelsdisplayed8( int m, int M)
        {

            this.odlm = m;
            this.odlM = M;
            double imagerange = Math.Abs(Convert.ToDouble(M - m));

            bmp = new Bitmap(Width,Height, System.Drawing.Imaging.PixelFormat.Format24bppRgb);

            BitmapData bmd = bmp.LockBits(new Rectangle(0, 0, Width,Height), System.Drawing.Imaging.ImageLockMode.ReadOnly, bmp.PixelFormat);

            // This 'unsafe' part of the code populates the bitmap bmp with data stored in pix8.
            // It does so using pointers, and therefore the need for 'unsafe'. 
            if (zeroiswhite == false)
            {
                unsafe
                {
                    int pixelSize = 3;
                    int i, j, j1;
                    byte bite;

                    for (i = 0; i < bmd.Height; ++i)
                    {
                        byte* row = (byte*)bmd.Scan0 + (i * bmd.Stride);

                        for (j = 0; j < bmd.Width; ++j)
                        {

                            bite = Convert.ToByte((Convert.ToDouble(pixels8[i * bmd.Width + j] - m) * 255.0 / imagerange));
                            if (bite > 255)
                                bite = 255;

                            if (bite < 0)
                                bite = 0;

                            j1 = j * pixelSize;
                            row[j1] = bite;            // Red
                            row[j1 + 1] = bite;            // Green
                            row[j1 + 2] = bite;            // Blue
                        }
                    }
                } // end unsafe      
                bmp.UnlockBits(bmd);

                //bmp.Dispose();
            }
            else 
            {
                unsafe
                {
                    int pixelSize = 3;
                    int i, j, j1;
                    byte bite;

                    for (i = 0; i < bmd.Height; ++i)
                    {
                        byte* row = (byte*)bmd.Scan0 + (i * bmd.Stride);

                        for (j = 0; j < bmd.Width; ++j)
                        {

                            bite = Convert.ToByte((Convert.ToDouble(M - pixels8[i * bmd.Width + j]) * 255.0 / imagerange));
                            if (bite > 255)
                                bite = 255;

                            if (bite < 0)
                                bite = 0;

                            j1 = j * pixelSize;
                            row[j1] = bite;            // Red
                            row[j1 + 1] = bite;            // Green
                            row[j1 + 2] = bite;            // Blue
                        }
                    }
                } // end unsafe      
                bmp.UnlockBits(bmd);

            }

        }

        /// <summary>
        /// Read the image pixels and create a bitmap for displaying the image.
        /// </summary>
        /// <param name="m">
        /// It's the smallest value of the desired value range to be displayed.
        /// </param>
        /// <param name="M">
        /// It's the biggest value of the desired value range to be displayed.
        /// </param>
        private  void changescaleofpixelsdisplayed16(int m, int M)
        {

            this.odlm = m;
            this.odlM = M;
            double imagerange = Math.Abs( Convert.ToDouble(M - m));

            bmp = new Bitmap(Width,Height, System.Drawing.Imaging.PixelFormat.Format24bppRgb);

            BitmapData bmd = bmp.LockBits(new Rectangle(0, 0, Width,Height), System.Drawing.Imaging.ImageLockMode.ReadOnly, bmp.PixelFormat);

            if (zeroiswhite == false)
            {
                unsafe
                {

                    double b;

                    for (int i = 0; i < bmd.Height; ++i)
                    {
                        byte* row = (byte*)bmd.Scan0 + (i * bmd.Stride);
                        for (int j = 0; j < bmd.Width; ++j)
                        {
                            b = Convert.ToDouble(pixels16[i * bmd.Width + j] - m) * 255.0 / imagerange;


                            if (b > 255)
                                b = 255;

                            if (b < 0)
                                b = 0;


                            row[j * 3] = (byte)b;            // Red
                            row[j * 3 + 1] = (byte)b;           // Green
                            row[j * 3 + 2] = (byte)b;           // Blue


                        }
                    }
                } // end unsafe   
                bmp.UnlockBits(bmd);
                //bmp.Dispose();
                //don't use this dispose I don't know why but it causes a strange error.
            }
            else 
            {
                unsafe
                {

                    double b;

                    for (int i = 0; i < bmd.Height; ++i)
                    {
                        byte* row = (byte*)bmd.Scan0 + (i * bmd.Stride);
                        for (int j = 0; j < bmd.Width; ++j)
                        {
                            b = Convert.ToDouble(M-pixels16[i * bmd.Width + j]) * 255.0 / imagerange;


                            if (b > 255)
                                b = 255;

                            if (b < 0)
                                b = 0;


                            row[j * 3] = (byte)b;            // Red
                            row[j * 3 + 1] = (byte)b;           // Green
                            row[j * 3 + 2] = (byte)b;           // Blue


                        }
                    }
                } // end unsafe   
                bmp.UnlockBits(bmd);
            }
        }

        /// <summary>
        /// Read the image pixels and create a bitmap for displaying the image.
        /// </summary>
        private void  changescaleofpixelsdisplayed8()
        {
            double imagerange = Math.Abs( Convert.ToDouble(max - min));

            bmp = new Bitmap(Width,Height, System.Drawing.Imaging.PixelFormat.Format24bppRgb);

            BitmapData bmd = bmp.LockBits(new Rectangle(0, 0, Width,Height), System.Drawing.Imaging.ImageLockMode.ReadOnly, bmp.PixelFormat);

            if (zeroiswhite == false)
            {
                unsafe
                {

                    double b;

                    for (int i = 0; i < bmd.Height; ++i)
                    {
                        byte* row = (byte*)bmd.Scan0 + (i * bmd.Stride);
                        for (int j = 0; j < bmd.Width; ++j)
                        {
                            b = Convert.ToDouble(pixels8[i * bmd.Width + j] - min) * 255.0 / imagerange;


                            if (b > 255)
                                b = 255;

                            if (b < 0)
                                b = 0;


                            row[j * 3] = (byte)b;            // Red
                            row[j * 3 + 1] = (byte)b;            // Green
                            row[j * 3 + 2] = (byte)b;            // Blue


                        }
                    }
                } // end unsafe   
                bmp.UnlockBits(bmd);
            }
            else
            {
                unsafe
                {

                    double b;

                    for (int i = 0; i < bmd.Height; ++i)
                    {
                        byte* row = (byte*)bmd.Scan0 + (i * bmd.Stride);
                        for (int j = 0; j < bmd.Width; ++j)
                        {
                            b = Convert.ToDouble(max - pixels8[i * bmd.Width + j] ) * 255.0 / imagerange;


                            if (b > 255)
                                b = 255;

                            if (b < 0)
                                b = 0;


                            row[j * 3] = (byte)b;            // Red
                            row[j * 3 + 1] = (byte)b;            // Green
                            row[j * 3 + 2] = (byte)b;            // Blue


                        }
                    }
                } // end unsafe   
                bmp.UnlockBits(bmd);
            }
        }

        /// <summary>
        /// Read the image pixels and create a bitmap for displaying the image.
        /// </summary>
        private void changescaleofpixelsdisplayed16()
        {
            double imagerange = Math.Abs( Convert.ToDouble(max - min));

            bmp = new Bitmap(Width, Height, System.Drawing.Imaging.PixelFormat.Format24bppRgb);

            BitmapData bmd = bmp.LockBits(new Rectangle(0, 0, Width, Height), System.Drawing.Imaging.ImageLockMode.ReadOnly, bmp.PixelFormat);

            if (zeroiswhite == false)
            {
                unsafe
                {

                    byte b;

                    for (int i = 0; i < bmd.Height; ++i)
                    {
                        byte* row = (byte*)bmd.Scan0 + (i * bmd.Stride);
                        for (int j = 0; j < bmd.Width; ++j)
                        {
                            b = (byte)(Convert.ToDouble(pixels16[i * bmd.Width + j] - min) * 255.0 / imagerange);

                            //not doing this saves about 30% of local ticks.
                            if (b > 255)
                                b = 255;

                            if (b < 0)
                                b = 0;


                            row[j * 3] = b;            // Red
                            row[j * 3 + 1] = b;            // Green
                            row[j * 3 + 2] = b;            // Blue

                            //Attention: You can insert your histogram here!
                        }
                    }
                } // end unsafe   
                bmp.UnlockBits(bmd);
            }
            else
            {
                unsafe
                {

                    byte b;

                    for (int i = 0; i < bmd.Height; ++i)
                    {
                        byte* row = (byte*)bmd.Scan0 + (i * bmd.Stride);
                        for (int j = 0; j < bmd.Width; ++j)
                        {
                            b = (byte)(Convert.ToDouble(max - pixels16[i * bmd.Width + j]) * 255.0 / imagerange);

                            //not doing this saves about 30% of local ticks.
                            if (b > 255)
                                b = 255;

                            if (b < 0)
                                b = 0;


                            row[j * 3] = b;            // Red
                            row[j * 3 + 1] = b;            // Green
                            row[j * 3 + 2] = b;            // Blue

                            //Attention: You can insert your histogram here!
                        }
                    }
                } // end unsafe   
                bmp.UnlockBits(bmd);
            }
        }


        /// <summary>
        /// Set the maximum and the minimum value for a 8 bit depth vector
        /// </summary>
        private void setmaxandmin8()
        {
            max = 0;
            min = 1000;

            // Compute imageMax and imageMin;

            for (int i = 0; i < pixels16.Length; i++)
            {
                if (pixels8[i] > max)
                    max = pixels8[i];

                if (pixels8[i] < min)
                    min = pixels8[i];
            }
        }


        /// <summary>
        /// Set the maximum and the minimum value for a 16 bit depth vector
        /// </summary>
        private void setmaxandmin16()
        {
            max = 0;
            min = 100000;

            for (int i = 0; i < pixels16.Length; i++)
            {
                if (pixels16[i] > max)
                    max = pixels16[i];

                if (pixels16[i] < min)
                    min = pixels16[i];

            }
        }

        /// <summary>
        /// Finds the max and min values inside a desired area and returns this value.
        /// </summary>
        /// <param name="calc"></param>
        /// <returns></returns>
        private Point maxandmin8(Rectangle calc)
        {
            int smax = 0;
            int smin = 100000;

            for (int j = 0; j < calc.Height; j++)
                for (int i = calc.X; i <= calc.X + calc.Width; i++)
                {
                    if (pixels8[i + (j + calc.Y) * Width] > smax)
                        smax = pixels8[i + (j + calc.Y) * Width];
                    if (pixels8[i + (j + calc.Y) * Width] < smin)
                        smin = pixels8[i + (j + calc.Y) * Width];
                }

            return new Point(smin, smax);

        }

        /// <summary>
        /// Finds the max and min values inside a desired area and returns this value.
        /// </summary>
        /// <param name="calc">
        /// Coordinates of the desired piece of image where the funcion will work.
        /// </param>
        /// <returns></returns>
        private Point maxandmin16(Rectangle calc)
        {
            int smax = 0;
            int smin = 100000;

            for (int j = 0; j < calc.Height; j++)
                for (int i = calc.X; i <= calc.X + calc.Width; i++)
                {
                    if (pixels16[i + (j + calc.Y) * Width] > smax)
                        smax = pixels16[i + (j + calc.Y) * Width];
                    if (pixels16[i + (j + calc.Y) * Width] < smin)
                        smin = pixels16[i + (j + calc.Y) * Width];
                }

            return new Point(smin, smax);

        }
        /*
        private Point maxandmin16overaline(Point p1, Point p2)
        {

            int smax = 0;
            int smin = 100000;

            //straight vertical line
            if (p1.X == p2.X)
            {
                if (p1.Y > p2.Y)
                {
                    for (int i = p2.Y; i < p1.Y; i++)
                        if (pixels16[p1.X + (i) * Width]> smax)
                            smax = pixels16[p1.X + (i) * Width];
                }
                else
                {
                    for (int i = p1.Y; i < p2.Y; i++)
                        if (pixels16[p1.X + (i) * Width]> smax)
                            smax = pixels16[p1.X + (i) * Width];
                }
            }
            if (p1.Y == p2.Y) 
            {

            }
        }

            /*
            {
                if (p2.Y < p1.Y)
                {
                    for (int i = 0; i < p2.X - p1.X; i++)
                    {
                        if (pixels16[i + p1.X + (int)(i*((float)(p1.Y-p2.Y)/(p2.X-p2.Y)) * Width)] > smax)
                            smax = pixels16[i + () * Width];
                        if (pixels16[i + () * Width] < smin)
                            smin = pixels16[i + () * Width];
                    }

                }
                else
                {

                }

            }
            else
            {
                if (p2.Y > p1.Y)
                {

                }
                else
                {

                }
            }

           

            return new Point(smin, smax);
        }
    */
        /// <summary>
        /// Finds the max and min values inside a desired area and returns this value.
        /// </summary>
        /// <param name="calc">
        /// Coordinates of the desired piece of image where the funcion will work.
        /// </param>
        /// <returns></returns>
        private double averagevalue8(Rectangle calc)
        {
            double sum = 0;

            for (int j = 0; j < calc.Height; j++)
                for (int i = calc.X; i < calc.X + calc.Width; i++)
                {
                    sum += pixels8[i + (j + calc.Y) * Width];
                }

            return sum / (calc.Width * calc.Height);
        }


       /// <summary>
        /// Calculates the average value of the pixels in the range of the vector.
       /// </summary>
       /// <param name="calc">
        /// Coordinates of the desired piece of image where the funcion will work.
        /// </param>
       /// <returns></returns>
        private double averagevalue16(Rectangle calc)
        {
            double sum = 0;

            for (int j = 0; j < calc.Height; j++)
                for (int i = calc.X; i < calc.X + calc.Width; i++)
                {
                    sum += pixels16[i + (j + calc.Y) * Width];
                }

            return sum / (calc.Width * calc.Height);
        }


        /// <summary>
        /// Calculates the stddev value of the pixels in the range of the vector.
        /// </summary>
        /// <param name="calc">
        /// Coordinates of the desired piece of image where the funcion will work.
        /// </param>
        /// <param name="average"></param>
        /// <returns></returns>
        private double stddev8(Rectangle calc, double average)
        {
            double std = 0;
            for (int j = 0; j < calc.Height; j++)
                for (int i = calc.X; i < calc.X + calc.Width; i++)
                {
                    std += Math.Pow((pixels8[i + (j + calc.Y) * Width] - average), 2);
                }

            return Math.Sqrt(std / (calc.Width * calc.Height - 1));
        }


        /// <summary>
        /// Calculates the stddev value of the pixels in the range of the vector.
        /// </summary>
        /// <param name="calc">
        /// Coordinates of the desired piece of image where the funcion will work.
        /// </param>
        /// <param name="average"></param>
        /// <returns></returns>
        private double stddev16(Rectangle calc, double average)
        {
            double std = 0;
            for (int j = 0; j < calc.Height; j++)
                for (int i = calc.X; i < calc.X + calc.Width; i++)
                {
                    std += Math.Pow((pixels16[i + (j + calc.Y) * Width] - average), 2);
                }

            return Math.Sqrt(std / (calc.Width * calc.Height - 1));
        }

        /// <summary>
        /// Calculates the differential unifformidity inside a specific area.
        /// </summary>
        /// <param name="calc"></param>
        private double unifordiff16(Rectangle calc)
        {

            int maxp, mimp;

            float udiff = 0;
            float udiffp;

            if (calc.Height >= 5)
                for (int j = 0; j < calc.Height - 3; j++)
                    for (int i = calc.X; i <= calc.X + calc.Width; i++)
                    {
                        maxp = 0;
                        mimp = 10000;

                        for (int k = 0; k < 5; k++)
                        {
                            if (pixels16[i + (j + calc.Y + k) * Width] > maxp)
                                maxp = pixels16[i + (j + calc.Y + k) * Width];
                            if (pixels16[i + (j + calc.Y + k) * Width] < mimp)
                                mimp = pixels16[i + (j + calc.Y + k) * Width];
                        }
                        udiffp = (float)(maxp - mimp) / (maxp + mimp);
                        if (udiffp > udiff)
                            udiff = udiffp;
                    }

            if (calc.Width >= 5)
                for (int j = 0; j < calc.Height; j++)
                    for (int i = calc.X; i <= calc.X + calc.Width - 3; i++)
                    {
                        maxp = 0;
                        mimp = 10000;

                        for (int k = 0; k < 5; k++)
                        {
                            if (pixels16[i + k + (j + calc.Y) * Width] > maxp)
                                maxp = pixels16[i + k + (j + calc.Y) * Width];
                            if (pixels16[i + k + (j + calc.Y) * Width] < mimp)
                                mimp = pixels16[i + k + (j + calc.Y) * Width];
                        }
                        udiffp = (float)(maxp - mimp) / (maxp + mimp);
                        if (udiffp > udiff)
                            udiff = udiffp;
                    }

            return  udiff *= 100;

        }


        /// <summary>
        /// Creates a smaller image based on the pixels that are shown at the screen.
        /// </summary>
        /// <param name="sender">
        /// Size is the desired squared width value.
        /// </param>
        /// <param name="e"></param>
        private void acumulateimage8(int size)
        {

            int mimifactor = this.pixels8.Length / size * size;
            ushort color;
            ushort[] newimagepixels = new ushort[size * size];


            for (int i = 0; i < pixels8.Length / mimifactor * mimifactor; i++)
                for (int j = 0; j < pixels8.Length / mimifactor * mimifactor; j++)
                {
                    color = 0;
                    for (int k1 = 0; k1 < mimifactor; k1++)
                        for (int k2 = 0; k2 < mimifactor; k2++)
                            color += pixels8[i * mimifactor + k1 + (j * mimifactor + k2) * Width];

                    newimagepixels[i + j * size] = color;
                }

            pixels16 = newimagepixels;
            setmaxandmin16();
            changescaleofpixelsdisplayed16(min, max);

        }


        /// <summary>
        /// Creates a smaller image based on the pixels that are shown at the screen.
        /// </summary>
        /// <param name="sender">
        /// Size is the desired squared width value.
        /// </param>
        /// <param name="e"></param>
        private void acumulateimage16( int size)
        {
            
            int mimifactor = this.pixels16.Length / size * size;
            ushort color;
            ushort[] newimagepixels = new ushort[size * size];


            for (int i = 0; i < pixels16.Length  / mimifactor*mimifactor; i++)
                for (int j = 0; j < pixels16.Length / mimifactor*mimifactor; j++)
                {
                    color = 0;
                    for (int k1 = 0; k1 < mimifactor; k1++)
                        for (int k2 = 0; k2 < mimifactor; k2++)
                            color += pixels16[i * mimifactor + k1+(j * mimifactor + k2) * Width];

                    newimagepixels[i + j * size] = color;
                }

            pixels16 = newimagepixels;
            setmaxandmin16();
            changescaleofpixelsdisplayed16(min, max);

        }


        /// <summary>
        /// A convolution kernel is applyed to the image (In the DICOM data).
        /// </summary>
        /// <param name="sender"></param>
        /// <param name="e"></param>
        private void filter16(int size)
        {
            ushort[] newpixels16 = new ushort[size * size];
            for (int i = 1; i < size - 1; i++)
                for (int j = 1; j < size - 1; j++)
                {
                    newpixels16[i + j * size] = pixels16[i + j * size];
                }



            for (int i = 1; i < size - 1; i++)
                for (int j = 1; j < size - 1; j++)
                {
                    if (pixels16[i + j * size] != 0)
                        newpixels16[i + j * size] =
                            (ushort)((float)
                            (pixels16[i - 1 + (j - 1) * size] + 2 * pixels16[i + (j - 1) * size] + pixels16[i + 1 + (j - 1) * size] +
                            2 * pixels16[i - 1 + (j) * size] + 4 * pixels16[i + (j) * size] + 2 * pixels16[i + 1 + (j) * size] +
                            pixels16[i - 1 + (j + 1) * size] + 2 * pixels16[i + (j + 1) * size] + pixels16[i + 1 + (j + 1) * size]) / 16);
                }

            pixels16 = newpixels16;

        }

        /// <summary>
        /// Calculates the average value in the center of the image.
        /// </summary>
        /// <param name="size">
        /// Size is the value of the square width.
        /// </param>
        /// <returns></returns>
        private float CFOVaverage(int size)
        {
            float average = 0;

            for (int i = size / 4; i < size / 2; i++)
                for (int j = size / 4; j < size / 2; j++)
                    average += pixels16[i + j * size];

            average = 16 * average / (size * size);

            return average;
        }

        /// <summary>
        /// It is a implementation of the NEMA standard. In this one the pixels sourrounding a zero pixel will be set to zero. 
        /// And the edges pixels with brightness less than 75% of the cfov average will be set to zero too.
        /// </summary>
        /// <param name="size">
        /// Size is the width of the squared image.
        /// </param>
        private void removebadpixels(int size)
        {

            ushort[] newpixels16 = new ushort[size * size];


            for (int i = 0; i < size; i++)
                for (int j = 0; j < size; j++)
                {
                    if (pixels16[i + j * size] == 0)
                    {
                        if (i - 1 > 0)
                            newpixels16[i - 1 + j * size] = 0;
                        if (i + 1 < size)
                            newpixels16[i + 1 + j * size] = 0;
                        if (j - 1 > 0)
                            newpixels16[i + (j - 1) * size] = 0;
                        if (j + 1 < size)
                            newpixels16[i + (j + 1) * size] = 0;
                    }
                    else
                        newpixels16[i + j * size] = pixels16[i + j * size];
                }

            pixels16 = newpixels16;


            float mean = CFOVaverage(size);

            for (int i = 0; i < size; i++)
                for (int j = 0; j < 15; j++)
                {
                    if (pixels16[i + j * size] < (ushort)(0.75 * mean))
                        pixels16[i + j * size] = 0;

                }

            for (int i = 0; i < size; i++)
                for (int j = size - 15; j < 64; j++)
                {
                    if (pixels16[i + j * size] < (ushort)(0.75 * mean))
                        pixels16[i + j * size] = 0;

                }

            for (int i = 0; i < 10; i++)
                for (int j = 0; j < size; j++)
                {
                    if (pixels16[i + j * size] < (ushort)(0.75 * mean))
                        pixels16[i + j * size] = 0;

                }
            for (int i = size - 10; i < size; i++)
                for (int j = 0; j < size; j++)
                {
                    if (pixels16[i + j * size] < (ushort)(0.75 * mean))
                        pixels16[i + j * size] = 0;

                }

        }


        /// <summary>
        /// It works as a image filter that inverts the brightness of each pixel.
        /// </summary>
        private void invertpixels8()
        {
            for (int i = 0; i < pixels8.Length; i++)
                pixels8[i] = (byte)(max - pixels8[i]);
        }

        /// <summary>
        /// It works as a image filter that inverts the brightness of each pixel.
        /// </summary>
        private void invertpixels16()
        {
            for (int i = 0; i < pixels16.Length; i++)
                pixels16[i] = (ushort)(max - pixels16[i]);
        }

        /// <summary>
        /// Develops a histogram equalization over the image pixels
        /// </summary>
        public void histogramequalizataion() 
        {
            if (bitDepth == 8)
            {
                int size = 2 << bitDepth;
                int[] h255 = new int[size];
                for (int i = 0; i < size; i++)
                    h255[i] = 0;

                //Basic Grayscale histogram
                for (int i = 0; i < pixels8.Length; i++)
                    h255[pixels8[i]]++;

                //Acumulated histogram
                for (int i = 0; i < size-1; i++)
                    h255[i + 1] += h255[i];

                for (int i = 0; i < pixels8.Length; i++)
                    pixels8[i] = (byte)(((h255[pixels8[i]] - h255[0]) * (size)) / (pixels8.Length - h255[0]));
            }

            if (bitDepth == 16)
            {
                int size = 4096;
                int[] h4095 = new int[size];
                for (int i = 0; i < size; i++)
                    h4095[i] = 0;

                //Basic Grayscale histogram
                for (int i = 0; i < pixels16.Length; i++)
                    h4095[pixels16[i]]++;

                //Accumulated histogram
                for (int i = 0; i < size - 1; i++)
                    h4095[i + 1] += h4095[i];


                //This section of the code is intended to find the maximum value and the minimum one over the cdf (accumulated histogram)
                int hmin = 1000000000, hmax = 0;

                for (int i = 0; i < size; i++)
                {
                    if (h4095[i] > hmax)
                        hmax = h4095[i];
                    if (h4095[i] < hmin)
                        hmin = h4095[i];
                }              

                //Change the value of the image pixels
                for (int i = 0; i < pixels16.Length; i++)
                    pixels16[i] = (ushort)((size - 1) * ((float)(h4095[pixels16[i]] - hmin) / (hmax - hmin)));
            }

        }

        public void convolution(int[,] kernel, Rectangle calc) 
        {

            if (kernel.GetLength(0) != kernel.GetLength(1))
                MessageBox.Show("O kernel selecionado não é quadrado!");
            else 
            {
                if (bitDepth == 16)
                    convolution16(kernel, calc);
                else
                    convolution8(kernel, calc);
            }

        }

        public dicomimage getimagerectangle(Rectangle calc)
        {
            ushort[] newrectangle = new ushort[calc.Width*calc.Height];
            dicomimage ndicom;
            if (bitDepth == 16)
            { 
                for (int i = 0; i < calc.Width; i++)
                    for (int j = 0; j < calc.Height; j++)
                    {
                        newrectangle[i+j*calc.Width] = this.pixels16[i+calc.X+ (j+calc.Y)*this.Width];
                    }
                ndicom = new dicomimage(newrectangle);
                ndicom.bitDepth = 16;
                ndicom.Width = calc.Width;
                ndicom.Height = calc.Height;

            }
            else 
            {
                for (int i = 0; i < calc.Width; i++)
                    for (int j = 0; j < calc.Height; j++)
                    {
                        newrectangle[i + j * calc.Width] = pixels8[i + calc.X + (j + calc.Y) * this.Width];
                    }
                ndicom = new dicomimage(newrectangle);
                ndicom.bitDepth = 16;
                ndicom.Width = calc.Width;
                ndicom.Height = calc.Height;
 
            }

            return ndicom;
 
        }

        public void convolution(int[,] kernel)
        {

            if (kernel.GetLength(0) != kernel.GetLength(1))
                MessageBox.Show("O kernel selecionado não é quadrado!");
            else
            {
                if (bitDepth == 16)
                    convolution16(kernel, new Rectangle(0, 0, this.Width, this.Height));
                else
                    convolution8(kernel, new Rectangle(0, 0, this.Width, this.Height));
            }

        }

        /// <summary>
        /// Creates a function that calculates the convolution between the image and a desired kernel.
        /// </summary>
        /// <param name="kernel"></param>
        /// <param name="calc"></param>
        private void convolution16(int[,] kernel, Rectangle calc)
        {
            ushort[] newpixels16 = new ushort[pixels16.Length];
            int partialresult;
            int normalization = 0;
            

            for (int i = 0; i < kernel.GetLength(0); i++)
                for (int j = 0; j < kernel.GetLength(1); j++)
                    normalization += kernel[i, j];

            if (normalization == 0)
                normalization = 1;


            for (int i = 0; i < calc.Width; i++)
                for (int j = 0; j < calc.Height; j++)
                {
                    partialresult = 0;
                    for (int k = -kernel.GetLength(0) / 2; k <= kernel.GetLength(0) / 2; k++)
                        for (int l = -kernel.GetLength(1) / 2; l <= kernel.GetLength(1) / 2; l++)
                            if ((i + k >= 0) && (i + k < Width) && (j + l >= 0) && (j + l < Height))
                                partialresult += pixels16[(i + k + calc.X) + (j + l + calc.Y) * Width] * kernel[k + kernel.GetLength(0) / 2, l + kernel.GetLength(1) / 2];
                    newpixels16[i + j * Width] = (ushort)Math.Abs((partialresult / normalization));
                }

            pixels16 = newpixels16;
        }

        /// <summary>
        /// Creates a function that calculates the convolution between the image and a desired kernel.
        /// </summary>
        /// <param name="kernel"></param>
        /// <param name="calc"></param>
        private void convolution8(int[,] kernel, Rectangle calc)
        {
            byte[] newpixels8 = new byte[pixels8.Length];
            int partialresult;
            int normalization = 0;


            for (int i = 0; i < kernel.GetLength(0); i++)
                for (int j = 0; j < kernel.GetLength(1); j++)
                    normalization += kernel[i, j];

            if (normalization == 0)
                normalization = 1;

            for (int i = 0; i < calc.Width; i++)
                for (int j = 0; j < calc.Height; j++)
                {
                    partialresult = 0;
                    for (int k = -kernel.GetLength(0) / 2; k <= kernel.GetLength(0) / 2; k++)
                        for (int l = -kernel.GetLength(1) / 2; l <= kernel.GetLength(1) / 2; l++)
                            if ((i + k >= 0) && (i + k < Width) && (j + l >= 0) && (j + l < Height))
                                partialresult += pixels8[(i + k + calc.X) + (j + l+calc.Y) * Width] * kernel[k + kernel.GetLength(0) / 2, l + kernel.GetLength(1) / 2];
                    newpixels8[i + j * Width] = (byte)Math.Abs((partialresult / normalization));
                }

            pixels8 = newpixels8;
        }

        /// <summary>
        /// This function creates a small image to be displayed (image preview).
        /// </summary>
        /// <param name="size">
        /// Sets the size (width and heigth) of the desired preview image.
        /// </param>
        /// <returns>
        /// Returns a pixel vector.
        /// </returns>
        public ushort[] preview16(int size) 
        {
            if (size < this.Width && size < this.Height)
            {
                int w = this.Width / size;
                int h = this.Height / size;
                ushort[] previewimage = new ushort[size * size];
                for (int i = 0; i < size; i++)
                    for (int j = 0; j < size; j++)
                        previewimage[i + j * size] = pixels16[i * w + j * Width * h];

                return previewimage;
            }
            else
                return null;
        }

        /// <summary>
        /// This function creates a small image to be displayed (image preview).
        /// </summary>
        /// <param name="size">
        /// Sets the size (width and heigth) of the desired preview image.
        /// </param>
        /// <returns>
        /// Returns a pixel vector.
        /// </returns>
        public byte[] preview8(int size)
        {
            if (size < this.Width && size < this.Height)
            {
                int w = this.Width / size;
                int h = this.Height / size;
                byte[] previewimage = new byte[size * size];
                for (int i = 0; i < size; i++)
                    for (int j = 0; j < size; j++)
                        previewimage[i + j * size] = pixels8[i * w + j * Width * h];

                return previewimage;
            }
            else
                return null;
        }


        public int[] acumulateoveravector(bool horizontal)
        {
            int[] line = new int[this.Width];

            if (bitDepth == 16)
            {
                if (horizontal)
                {
                    for (int i = 0; i < this.Width; i++)
                        for (int j = 0; j < this.Height; j++)
                        {
                            line[i] += pixels16[i +  + (j ) * this.Width];
                        }
                }
                else
                {
                    line = new int[this.Height];
                    for (int i = 0; i < this.Height; i++)
                        for (int j = 0; j < this.Width; j++)
                        {
                            line[i] += pixels16[j  + (i ) * this.Width];
                        }

                }
            }
            else 
            {
                if (horizontal)
                {
                    for (int i = 0; i < this.Width; i++)
                        for (int j = 0; j < this.Height; j++)
                        {
                            line[i] += pixels8[i  + (j) * this.Width];
                        }
                }
                else
                {
                    line = new int[this.Height];
                    for (int i = 0; i < this.Height; i++)
                        for (int j = 0; j < this.Width; j++)
                        {
                            line[i] += pixels8[j + (i) * this.Width];
                        }

                }

            }
            return line;
        }

        public void changecolortype(bool type) 
        {
            if (type == true)
            {
                if (pixels16 != null)
                    for (int i = 0; i < this.Width * this.Height; i++)
                        pixels16[i] = (ushort)(this.max - pixels16[i]);
                else
                    for (int i = 0; i < this.Width * this.Height; i++)
                        pixels8[i] = (byte)(this.max - pixels8[i]);
            }
            else
            {
                if (pixels16 != null)
                    for (int i = 0; i < this.Width * this.Height; i++)
                        pixels16[i] = (ushort)(this.max - pixels16[i]);
                else
                    for (int i = 0; i < this.Width * this.Height; i++)
                        pixels8[i] = (byte)(this.max - pixels8[i]);
 
            }
        }


        public double[] returnbrightprofilevertical(Rectangle start) 
        {
            if (this.bitDepth == 16)
                return returnbrightprofilevertical16(start);
            else
                return returnbrightprofilevertical8(start); 

        }

        private double[] returnbrightprofilevertical16(Rectangle start) 
        {
            int average;
            double[] vector = new double[start.Height];

            for (int j = 0; j < start.Height; j++)
            {
                average = 0;
                for (int i = start.X; i < start.X + start.Width; i++)
                {
                    average += pixels16[i + (j + start.Y) * Width];
                }
                vector[j] = average / start.Width;
            }

            return vector;
        }


        private double[] returnbrightprofilevertical8(Rectangle start)
        {
            int average;
            double[] vector = new double[start.Height];

            for (int j = 0; j < start.Height; j++)
            {
                average = 0;
                for (int i = start.X; i < start.X + start.Width; i++)
                {
                    average += pixels8[i + (j + start.Y) * Width];
                }
                vector[j] = average / start.Width;
            }

            return vector;
        }

        public double[] returnbrightprofilehorizontal(Rectangle start)
        {
            if (this.bitDepth == 16)
                return returnbrightprofilehorizontal16(start);
            else
                return returnbrightprofilehorizontal8(start);

        }

        private double[] returnbrightprofilehorizontal16(Rectangle start)
        {
            int average;
            double[] vector = new double[start.Width];

            for (int i = start.X; i < start.X + start.Width; i++)
            {
                average = 0;
                for (int j = 0; j < start.Height; j++)
                {
                    average += pixels16[i + (j + start.Y) * Width];
                }
                vector[i - start.X] = average / start.Height;
            }

            return vector;
        }


        private double[] returnbrightprofilehorizontal8(Rectangle start)
        {
            int average;
            double[] vector = new double[start.Height];

            for (int i = start.X; i < start.X + start.Width; i++)
            {
                average = 0;
                for (int j = 0; j < start.Height; j++)
                {
                    average += pixels8[i + (j + start.Y) * Width];
                }
                vector[i - start.X] = average / start.Height;
            }

            return vector;
        }

    }
}
